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ABSTRACT 
In  this  thesis  results  of  a  computational  and  experimental 
investigation  of  the  aerodynamic  characteristics  of  a  sail 
section  used  in  windsurfing  sails  are  presented.   State-of-the- 
art  computational  methods  (panel,  direct  boundary  layer,  viscous - 
inviscid  interaction,  Euler,  and  steady/unsteady  Navier- Stokes) 
were  used  to  predict  the  aerodynamic  loading  and  stall 
characteristics.   These  predictions  were  found  to  be  in 
satisfactory  agreement  with  tuft  and  smoke  flow  visualization 
experiments  carried  out  in  the  Naval  Postgraduate  School  low 
speed  wind  tunnel  at  a  Reynolds  number  of  800,000.   Further,  all 
computational  work  was  completed  on  the  Silicon  Graphics  Indigo 
workstation  to  demonstrate  that  only  modest  computer  facilities 
will  be  necessary  for  these  methods  to  migrate  to  the  field  of 
sail  design. 
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I .   INTRODUCTION 

The  field  of  computational  fluid  dynamics  (CFD)  has  come  of 
age  over  the  past  several  years.  Rapid  advances  in  high  speed 
digital  computers,  their  subsequent  lower  cost  and  the  maturation 
of  CFD  codes  suggest  that  CFD  methods  should  spread  to  areas  beyond 
the  aerospace  industry  in  the  near  future.  The  discipline  of 
competitive  sailing  is  a  promising  area  that  could  benefit  from  the 
use  of  CFD  techniques.  Some  have  suggested  that  such  applications 
would  be  trivial.  However,  experience  shows  that  sailing  is  a 
fiercely  competitive  sport  ranging  from  amateur  levels  to  Olympic, 
professional  and  the  America's  Cup  levels.  In  practical  terms,  the 
bottom  line  has  always  been  performance,  which  directly  relates  to 
the  incorporation  of  technology  into  the  design  process.  The  aim 
of  this  research  will  be  to  demonstrate  that  CFD  techniques  may  be 
readily  adapted  to  the  field  of  sail  design. 

From  an  engineering  standpoint,  sailing  craft  are  challenging 
due  to  the  interaction  or  coupling  of  two  different  fluid  media, 
the  water's  surface  and  the  atmosphere.  Additionally,  the  two 
media  are  subject  to  separate  noise  or  turbulence  in  the  form  of 
waves,  wind  gradient  and  gusts.  The  objective  of  this  thesis  will 
be  to  apply  CFD  techniques  to  the  analysis  of  the  above  water 
portion  of  the  problem  and  to  provide  some  experimental 
verification.  The  scope  of  the  analysis  will  be  limited  to  a  two- 
dimensional  approach.  A  full  three-dimensional  analysis  would  have 
been  optimum,  but  the  complexity  of  the  problem  and  computing  power 
required  would  have  become  excessive. 


There  are  numerous  rig  and  sail  configurations  which  ought  to 
be  analyzed.  Planing  hull  sailing  craft,  in  general,  are  more 
challenging  because  they  are  capable  of  higher  speeds  than 
displacement  hull  yachts.  Planing  hull  craft,  moreover,  experience 
a  pronounced  rise  in  drag  prior  to  transitioning  from  a 
displacement  hull  mode  to  a  planing  hull. 

For  this  study  a  sailboard  or  windsurfer  rig  was  selected  over 
other  rig  and  sail  combinations  for  several  reasons.  The  primary 
reason  for  the  selection  of  this  rig  and  sail  combination  is  the 
rig's  simplicity,  Figure  1.1.  Its  shape  is  approximate  to  a  simple 
wing  or  foil  shape  of  moderate  aspect  ratio  and  is  favorable  when 
compared  to  a  sloop  or  other  multi-sail  configuration.  The 
acceleration  through  the  above  mentioned  planing  transition  is  of 
high  interest  for  study.  This  acceleration  or  quick  planing 
ability  is  of  paramount  importance  to  the  competitive  sailor  on  any 
planing  hull  dingy  or  sailboard.  To  accelerate  quickly  the  sailor 
in  effect  must  generate  large  lift  coefficients  from  the  sail. 
This  is  achieved  through  a  large  sheeting  angle  (high  angle  of 
attack)  in  combination  with  the  highly  cambered  shape  built  into 
the  sail  section  and  in  the  case  of  the  board  sailor  through  the 
dynamic  sheeting  or  pumping  of  the  sail.  The  effects  of  pumping 
are  well  known  to  all  competitive  sailors.  However,  this  sailing 
technique  is  strictly  prohibited  in  all  classes  except  some  types 
of  boardsailing.  The  pumping  method  utilized  is  a  refined  art  to 
many  board  sailors  to  the  degree  that  sail  designers  and 
competitors  note  this  quality  in  sails. 


Figure  1.1  Windsurfing  Sail 


The  second  regime  in  which  the  sail  performance  is  critical 
occurs  after  the  craft  has  fully  accelerated  on  a  plane.  As  a 
craft  accelerates  the  apparent  wind  shifts  forward  from  the 
sailor's  frame  of  reference  while  the  true  wind  velocity  remains 
constant,  Figure  1.2.  While  sailing  with  the  apparent  wind  far 
forward  the  sail  is  trimmed  at  or  near  the  maximum  lift  to  drag 
ratio.  The  maximum  lift  to  drag  ratio  in  fact  limits  how  close  to 
the  wind  a  vessel  may  sail  or  the  minimum  apparent  wind  angle  that 
may  be  sailed. 

The  technical  challenge  for  analysis  in  this  regime  of 
aerodynamics  lies  with  the  low  Reynolds  number  (Re)  below  one 
million.  The  bulk  of  data  to  date  for  airfoil  sections  is  for 
Reynolds  numbers  above  one  million.   CFD  work  has  been  done  at 
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Figure  1.2   Apparent  Wind  and  L/D  Relation 


Reynolds  number  below  one  million,  but  the  work  has  been  frequently 
restricted  to  internal  flows  in  turbines  or  ducts  or  to  research 
for  high  altitude  aircraft  where  the  Mach  numbers  are  significant. 
The  low  Reynolds  number  implies  that  the  viscous  forces  will  be 
large  and  require  thorough  investigation.  This  will  require  the 
use  of  viscid  codes  in  addition  to  simplified  inviscid  codes  to 
understand  the  nature  of  the  flow  field.  On  the  positive  side,  the 
problem  is  somewhat  simplified  in  that  it  is  entirely 
incompressible.  Physical  limitations  on  the  materials  used  and  the 
method  in  which  sails  are  constructed  produce  irregularities  and 
rough  surfaces  near  the  leading  edge  of  the  sections.  These 
imperfections,  while  not  desirable,  allow  for  the  assumption  that 
the  flow  is   fully  turbulent.   This  assumption  precludes  the 


requirement  for  and  validation  of  a  boundary  layer  transition 
model . 


II.     COMPUTATIONAL    APPROACH 

A.  GENERAL 

Computational  fluid  dynamics  has  been  largely  restricted  to 
well  funded  aerospace  research  and  development  groups  due  to  the 
need  for  access  to  computers  with  sufficient  power,  speed  and 
graphics  capability.  However,  recent  advances  in  computer  hardware 
have  eliminated  the  requirement  for  access  to  a  super  computer  to 
conduct  CFD  studies.  To  demonstrate  that  CFD  technology  is  mature 
enough  to  migrate  to  the  field  of  sail  design,  the  computational 
work  for  this  study  was  conducted  on  Silicon  Graphics  Indigo  and 
Iris  workstations  utilizing  Fortran  codes  and  Plot  3D  for  graphics 
presentation.  This  generation  of  workstations  have  largely 
redefined  the  conventional  boundaries  between  workstations, 
mainframes  and  mini-supercomputer  systems.  The  level  of  computing 
below  the  Unix  workstation  has  also  progressed  rapidly  with  the 
latest  generation  of  personal  computers  (PC)  built  around  the  Intel 
80486  processors.  These  computers  probably  have  sufficient  memory 
and  speed  for  this  application.  However,  at  this  time  they  are 
limited  by  a  lack  of  software  which  has  been  widely  available  for 
the  Unix  operating  system.  It  may  be  possible  in  the  near  future 
that  new  operating  systems  for  the  PC  will  include  the  necessary 
software  features.  The  Unix  based  workstation,  while  a  step  above 
the  PC  is  still  well,  within  the  resources  of  the  sail  design  and 
manufacture  business. 

B.  SAIL  SHAPE  GENERATION  TECHNIQUE 

To  commence  the  study  of  an  airfoil  the  first  requirement 
encountered  is  a  precise  definition  of  the  airfoil  section.  To 
proceed  with  the  analysis  of  a  sail  section  it  is  necessary  to 
develop  a  systematic  method  of  describing  sail  shapes  with  a 
similar  high  degree  of  precision. 

Sail  design  has  never  been  an  exact  science.  However,  sail 
designers  and  sailors  have  always  given  attention  to  the  same 
parameters  that  we  see  in  airfoil  design.   The  maximum  camber  and 


its  location  aft  of  the  leading  edge  are  the  two  primary  parameters 
used  to  describe  a  sail's  shape  (the  camber  of  a  sail  is  most 
commonly  referred  to  as  the  'draft'  by  the  sailing  community) . 
Several  additional  parameters  are  used  to  further  describe  the 
windsurfing  sail  section.  The  mast  (vertical  load  bearing  spar) 
has  a  cylindrical  section  internal  to  the  leading  edge  or  luff 
portion  of  the  sail.  The  mast  radius  is  directly  analogous  to  the 
leading  edge  radius  used  in  the  definition  of  many  airfoil 
sections.  From  the  mast  section  aft  two  flat  panels  are  used  to 
fair  the  mast  cross  section  to  the  thin  segment  of  the  section  that 
compose  the  majority  of  the  sail.  The  length  of  this  faring  has 
been  defined  as  the  luff  pocket  length. 

A  fortran  code,  'shape. f  Appendix  A,  was  developed  to 
precisely  define  a  sail  section  based  on  the  geometric  parameters 
mentioned  above.  Two  models  for  the  camber  line  generation  have 
been  included  in  the  program.  For  a  large  number  of  sails  the 
camber  line  can  be  approximated  by  a  circular  arc  aft  of  the 
maximum  camber  location  and  with  a  second  order  polynomial  forward. 
This  model  is  the  first  camber  option  for  sail  generation  which 
requires  the   following  input  arguments: 

mast  radius     (2.0%  nominal  value) 

luff  pocket  length 

thickness 

maximum  camber  or  draft 

maximum  camber  location 
To  describe  a  more  unique  section  with  flat  segments  or 
reflexed  areas  a  higher  order  polynomial  fit  has  been  included  as 
the  second  option  for  camber  generation.  For  this  routine  the 
camber  line  is  defined  by  points  through  which  a  Lagrange 
polynomial  is  fitted.  The  first  three  input  parameters  for  this 
option  remain  the  same  but  in  place  of  the  maximum  camber  and 
location  the  user  enters  n  number  of  x,y  coordinates  to  define  the 
camber  shape. 

'Shape. f  builds  the  section  from  the  common  trailing  edge 
point  working  forward.  The  generation  routine  uses  the  selected 


camber  model  to  describe  the  sail  camber  shape  aft  of  the  luff 
dimension  with  an  increment  to  add  the  prescribed  thickness.  From 
the  luff  dimension  forward  a  straight  line  is  placed  tangent  tHxthe 
mast  radius  to  describe  the  shape.  Seventy  panels  top  and  bottom 
from  the  trailing  edge  point  are  used  to  this  point  and  a  final  ten 
panels  are  used  to  depict  the  mast  radius.  The  151  x,y  coordinates 
that  form  the  basic  geometry  of  the  section  are  written  into  a  file 
'shape. out. ' 

A  final  segment  of  'shape. f  adds  thirty  wake  points  extended 
tangent  to  the  trailing  edge  for  grid  generation  use.  The  length 
between  the  wake  points  begins  with  the  same  initial  dx  distance 
with  each  subsequent  length  expanded  by  a  factor  of  1.1.  This 
second  set  of  data  points  including  the  wake  points  is  written  to 
'hypgen.out'  for  subsequent  grid  generation  use. 

The  primary  shape  used  for  the  CFD  evaluation  and  tunnel 
testing  is  depicted  in  Figure  2.1.  This  section  shape  was  typical 
of  the  first  fully  battened  sails  appearing  in  the  mid-1980's.  The 
shapes  used  at  that  time  were  generally  highly  cambered  with  the 
maximum  camber  well  forward  of  the  mid-chord  point.  The  shape 
depicted  is  irregular,  in  particular  the  forward  segment,  when 
compared  to  a  conventional  airfoil.  Obviously,  a  better  shape  from 
an  aerodynamic  perspective  would  be  desirable  but  physical 
limitations  with  the  mast,  materials  and  manufacture  prevent  the 
building  of  'ideal'  sections. 

Sail  shapes  have  slowly  evolved  through  a  process  of  trial  and 
error  into  shapes  similar  to  that  in  Figure  2.2.  The  current  shapes 
have  less  camber  with  the  maximum  camber  point  closer  to  the  mid- 
cord.  The  shapes  shown  in  Figures  2.1/2  were  both  generated  using 
the  circular  arc  and  second  order  polynomial  method.  The 
'shape. out'  files  containing  the  x,y  stations  for  the  two  sections 
are  contained  in  Appendix  A. 

Current  interest  among  sail  makers  has  been  to  make  the  aft 
sections  increasingly  flat.  The  shapes  may  promise  to  perform 
better.  However,  the  flat  areas  are  difficult  to  build  into  the 
sail  and  are  prone  to  flutter.  An  advanced  shape  similar  to  the  ones 


Figure  2.1   11.5%  Camber  Section 


Figure  2.2   8.5%  Camber  Section 
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currently  being  tested  was  derived  utilizing  the  Lagrange 
polynomial  generation  method  and  is  shown  in  Figure  2.3. 

C.  PANEL  CODE 

Panel  codes  represent  an  introductory  means  to  model  the  flow 
field  around  a  body.  They  are  the  simplest  and  easiest  method  to 
explore  the  flow  field  for  the  subject  of  this  thesis.  Of  note, 
this  family  of  codes  can  easily  run  solutions  and  be  plotted  on  a 
modest  personal  computer. 

To  explain  the  panel  method  we  first  examine  Laplace's 
equation  governing  incompressible,  irrotational,  and  invisid  fluid 
flow: 

*  +*  =o  (2-1) 

xx         yy 

This  expression  is  a  homogeneous  linear  second  order  partial 
differential  equation.  The  linearity  of  the  equation  allows  the 
use  of  the  superposition  principle  to  describe  the  flow  field  with 
elementary  flow  elements.  For  the  case  of  a  two  dimensional 
airfoil  in  a  uniform  flow  the  system  may  be  expressed  as  a 
combination  of  uniform,  source  and  vortex  potential  elements: 

where 

4>.=Y-  (xcosa+ysina) 

(2.3) 


*„=A 
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Figure 


Advanced  Section 
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If  the  airfoil  section  is  decomposed  into  a  set  of  n  panels  to 
describe  the  section,  the  flow  field  can  be  represented  by  the 
superposition  of  n  sources  located  at  the  mid-panel  points,  a 
uniform  flow  and  a  constant  vortex  strength.  Hence,  the  total 
number  of  unknowns  is  n+1. 

By  applying  boundary  conditions  a  solution  to  the  system  may 
be  obtained.  The  flow  tangency  condition  at  the  panel  mid-points 
results  in  n  equations.  A  final  boundary  condition  is  arrived  at 
by  applying  the  Kutta  condition  forcing  the  velocities  of  the  upper 
and  lower  trailing  edge  panels  to  be  equal.  The  system  of 
elementary  potential  flows  is  now  reduced  to  n+1  equations  and 
unknowns  which  can  be  expressed  in  matrix  form: 

[A][q]=[B] 
where,  A  is  the  influence  coefficient  matrix,  q  is  a  column  vector 
with  the  values  of  the  n  source  strengths  and  the  vortex  strength. 
B  is  also  a  column  vector  equating  the  angular  difference  between 
the  free  stream  angle  of  attack  and  the  mid-point  tangents  and  the 
wake  condition. 

The  coding  for  a  panel  code  is  fairly  straightforward  but 
great  care  must  be  utilized  in  the  geometry  conventions  used  to 
specify  the  influence  coefficient  matrix.  The  code  'panel. f 
contained  in  Appendix  A  was  written  by  the  author  and  adapted  for 
use  on  sail  sections.  To  adequately  describe  the  irregular 
sections  of  interest  a  high  number  of  panels  was  required  with  the 
upper  limit  currently  200.  The  airfoil  shape  is  entered  as  an 
input  file  in  one  of  two  formats  by  the  user.  The  input  file  can 
consist  of  standard  x,y  ordinates  or  be  in  the  form  of  a  two 
dimensional  Plot-3D  grid  file.  If  a  Plot-3d  file  is  used  the 
upper  and  lower  trailing  edge  points  along  the  i  direction  must 
also  be  provided. 

From  a  computational  standpoint,  the  main  weakness  in  solving 
Laplace's  equation  occurs  in  the  'panel. f  code  when  inverting  the 
A  matrix.  The  thickness  of  the  section  and  number  of  panels 
utilized  have  a  strong  effect  on  the  error  when  inverting  the 
influence  coefficient  matrix.   The  actual  thickness  of  the  sail 
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section  is  of  the  order  of  five  mil's,  a  very  small  thickness, 
which  is  not  realistic  for  the  code's  precision.  Experience  has 
shown  0.5%  thickness  to  be  a  reasonable  compromise.  The  complex 
geometries  also  suggest  that  a  large  number  of  point  should  be  used 
to  describe  the  section  with  151  panels  having  been  determined  to 
be  sufficient. 

D.  BOUNDARY  LAYER  CODES 

The  panel  code  provided  the  simplest  method  to  obtain  results 
for  the  pressure  coefficient.  Similarly,  it  is  highly  desirable  to 
obtain  viscous  results  in  a  simple  and  timely  manner.  Boundary 
layer  codes  represent  the  next  logical  step  toward  understanding 
the  viscous  behavior  for  the  low  Reynolds  number  present. 

Two  boundary  layer  codes  were  explored  to  investigate  their 
suitability  for  sail  analysis.  The  first  of  the  codes,  'dbl2.f  a 
direct  boundary  layer  code,  was  employed  with  very  limited  utility. 
The  code  employs  a  panel  method  routine  to  compute  the  pressure 
distribution  which  is  then  used  to  compute  the  boundary  layer 
profiles.  The  second  code  utilized  was  the  viscous  invisid 
interaction  code  developed  by  Cebeci  at  McDonnell-Douglas  Aircraft 
Company.  This  code  carries  the  process  used  in  the  direct  boundary 
layer  code  a  stage  further.  After  computing  the  boundary  layer  the 
code  then  adds  the  additional  thickness  of  the  viscous  layer  to  the 
airfoil  shape.  This  new  effective  shape  is  then  run  through  the 
inviscid  scheme  again  to  compute  a  new  boundary  layer.  This 
iterative  routine  is  repeated  until  the  change  in  the  boundary 
layer  becomes  sufficiently  small.  The  iterative  boundary  layer 
technique  was  nothing  short  of  a  total  failure.  This  was  due  to 
the  inability  of  the  Smith-Hess  panel  routine  used  by  the  code  to 
successfully  handle  the  irregular  sections.  However,  some  insight 
into  the  flow  separation  characteristics  was  gained  through  this 
code. 

It  was  hoped  that  these  two  codes  or  methods  would  be  of  great 
value  as  they  represent  the  next  order  of  sophistication  above  the 
simple  potential  inviscid  solution.   Both  schemes  proved  to  be  a 
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disappointment  toward  the  analysis  of  the  sail  shapes.  While  these 
codes  work  well  for  standardized  airfoil  shapes  they  proved  to  be 
of  little  value  for  the  irregular  shapes.  Their  failure  is  due  to 
the  fact  that  they  do  not  handle  separated  flow  regions  well. 
Potential  flow  modeling  of  the  highly  irregular  sail  shapes 
produced  large  pressure  perturbations  which  caused  the  boundary 
layer  to  separate  in  the  numerical  solution.  After  several 
attempts  to  improve  on  the  solutions  these  methods  were  abandoned 
in  favor  of  the  Navier-Stokes  and  Euler  methods. 

E.  EULER  and  NAVIER-STOKES  METHOD 

The  Navier-Stokes  (NS)  equations  represent  the  most  robust 
tool  currently  in  use  in  the  field  of  computational  fluid  dynamics. 
The  derivation  of  the  NS  equations  and  their  CFD  solution  method 
will  not  be  discussed  because  they  are  well  documented  in 
References  1  and  5.  The  two-dimensional  vector  form  of  the 
equations  may  be  expressed  as: 


where, 


dt     dx      dy 


Q=(p,  pu,  pv,  e)  T 
pu 


(2.4) 
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The  Euler  solution  is  readily  obtained  from  the  NS  equations  when 
the  viscous  terms  are  discarded. 

The  Navier-Stokes  (NS)  code  used  to  examine  the  flow 
characteristics  of  the  sail  sections,  'ns2.f  was  developed  by 
Professor  J.  A.  Ekaterinaris  of  the  Navy-NASA  Joint  Institute  of 
Aeronautics.  Slight  modifications  to  'ns2.f  were  made  in  the  form 
of  additional  write  statements  to  save  unsteady  motion  solutions  at 
regular  time  intervals  and  to  simplify  steady  solution  inputs.  A 
call  to  the  grid  rotation  subroutine  was  added  to  preclude  the 
input  of  a  rotated  grid  when  running  'ns2.f.  The  major  features 
of  the  code  are: 

*  Upwind  differencing 

*  Baldwin-Lomax  turbulence  modeling 

*  Ability  to  restart  the  code 

Navier-Stokes  solutions  were  obtained  after  approximately  four 
thousand  iterations  from  a  uniform  flow  field,  based  on  the  density 
residuals,  Figure  2.4.  Solutions  could  be  obtained  with  a  fraction 
of  the  four  thousand  iterations  if  a  restart  was  initiated  from  a 
previous  solution  (example:  an  8  degree  angle  of  attack  solution  is 
used  to  compute  a  ten  degree  solution)  .  The  NS  solutions  were 
computed  using  the  Iris  Indigo  work  stations  in  the  CFD  laboratory. 
This  required  approximately  twelve  hours  of  CPU  time  for  four 
thousand  iterations  using  a  251  by  71  grid  run  with  a  Courant 
number  of  2100.  While  this  may  seem  to  be  a  huge  amount  of 
computer  time,  it  must  be  kept  in  mind  that  this  was  only  a  small 
work  station  and  that  few  of  the  solutions  required  four  thousand 
iterations  by  using  the  restart  feature.  Euler  solutions  can  also 
be  obtained  from  this  code  in  a  similar  manner  with  the  correct 
switches  in  the  input  file  to  discard  the  viscous  terms. 
F.   GRID  GENERATION 

The  need  for  a  grid  system  to  define  a  flow  field  around  a 
body  arises  from  the  necessity  of  transforming  from  the  physical 
domain  to  a  discrete  computational  domain.  The  theory  of  grid 
generation  will  not  be  covered  in  this  thesis.  Merkle's  text, 
Reference  5,  explains  this  process  as  do  several  other  texts  in  the 
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field  of  CFD.  Grid  generation  for  the  section  depicted  in  Figures 
2.1  and  2.2  was  accomplished  utilizing  the  code  'hypgen'. 

The  grid  type  selected  for  this  analysis  was  the  'C  grid 
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Figure  2.4  Navier-Stokes  Solution  Convergence 
which  is  most  commonly  used  for  2-dimensional  airfoils.  The  'C  is 
in  reference  to  the  shape  or  manner  in  which  this  type  of  grid  is 
wrapped  around  the  airfoil.  The  i  stations  are  along  the  wake  and 
airfoil  directions  while  the  k  stations  extend  from  the  airfoil 
surface  to  the  far  field.  The  primary  difficulty  encountered  was 
ensuring  the  orthogonality  of  the  grid  lines  around  the  leading 
edge  and  forward  portion  of  the  sail  sections.'  Additional  points 
along  the  airfoil  surface  were  added  to  the  luff  section  to  smooth 
the  interval  distance  along  the  i  direction  near  the  leading  edge. 
For  the  details  in  the  grid  generation  software  the  'hypgen'  user 
manual  should  be  consulted. 

The  full  dimensions  of  the  grids  used  are  shown  in  Figure  2.5. 
These  dimensions,  10  chords  lengths  ahead,  above,  below  and  30 
chord  lengths  aft  of  the  sail  section,  were  selected  to  ensure  that 
solutions  would  smoothly  match  the  far  field  boundary  conditions. 
The  grid  used  for  the  inviscid  Euler  solutions  is  shown,  Figure 
2.6.  For  the  NS  calculations  the  number  of  stations  from  the 
airfoil  to  the  far  field  was  increased  from  41  to  71  station.  This 
was  done  to  accurately  resolve  the  boundary  layer  near  the  airfoil 
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Figure  2.5   Full  Grid  Plot 
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surface.   The  increased  resolution  can  be  seen  in  Figure  2.7 


19 


Figure  2.6   Inviscid  Euler  Grid 
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Figure  2.7      Navier-Stokes   Grid 
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III.       EXPERIMENTAL    APPROACH 

A.  BACKGROUND 

The  most  frequent  criticism  of  computational  fluid  dynamic 
solutions  is  that  they  may  not  accurately  reflect  the  actual  flow 
present.  It  was  therefore  highly  desirable  to  validate  the 
computational  work  conducted  with  experimental  data  for  the 
sections  of  interest.  Investigations  into  previous  work  in  this 
field  failed  to  find  an  adequate  description  of  the  flow  around 
windsurfing  sail  sections.  The  majority  of  the  research  conducted 
in  this  field  has  been  for  sail  and  rig  combinations  for  specific 
yacht  types.  The  wind  tunnel  experiments  for  this  thesis  were 
designed  with  the  emphasis  on  predicting  the  separation  regions 
present  since  this  is  the  major  challenge  for  viscous  CFD  methods. 

B.  WIND  TUNNEL  EXPERIMENTAL  EQUIPMENT 

The  tunnel  selected  for  carrying  out  the  experiment  was  the  32 
by  45  inch  low  speed  wind  tunnel  located  at  NPS  in  Halligan  Hall. 
The  primary  reason  for  the  selection  of  this  tunnel  was  that  it  is 
capable  of  sufficient  velocities  to  run  tests  at  the  actual 
Reynolds  numbers  present  in  the  sailing  environment.  The  tunnel  is 
also  sufficiently  large  to  allow  for  reasonably  sized  models.  It 
was  desired  to  mount  the  model  sections  vertically  to  take 
advantage  of  the  large  optical  windows  on  the  sides  of  the  test 
sections  for  viewing  tufts.  In  addition,  by  using  the  greater 
tunnel  dimension  perpendicular  to  the  model  rotation  axis,  the 
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blockage  for  the  experiments  was  held  to  under  twelve  percent  for 
a  model  of  eighteen  inch  chord  with  up  to  eighteen  degrees  angle  of 
attack.  Reference  8  contains  a  detailed  description  of  the  tunnel 
and  its  use. 

The  next  area  addressed  was  the  choice  of  materials  to 
construct  the  models  from.  Particular  attention  was  given  to  the 
fact  that  the  airfoils  to  be  tested  have  very  thin  sections  over 
the  majority  of  the  chord.  A  conventional  rigid  model  would  have 
had  to  be  constructed  of  metal  or  fiberglass/composite  materials 
for  the  aft  section  and  of  wood  for  the  mast  and  luff  portions. 
This  approach  was  considered  but  abandoned  due  to  the  complexity, 
cost  and  long  lead  times  in  building  such  models. 

The  other  option  was  to  assemble  the  model  from  the  same 
materials  as  used  on  the  actual  sails.  This  approach  simplified 
the  model  building  process  and  allowed  for  easy  modification  of  the 
models.  The  use  of  the  same  materials  and  construction  techniques 
also  precluded  trying  to  simulate  the  roughness  associated  with  the 
different  fabrics,  films  and  seams  that  exist  on  windsurfing  sails. 

After  deciding  to  use  actual  2-dimensional  sails  as  models,  a 
means  to  hold  the  sail  and  mount  the  sail  in  the  test  section  was 
developed.  A  rig  was  designed  and  built  of  aluminum,  as  is 
depicted  in  Figure  3.1.  This  rig  allows  for  the  model  to  be  in 
tension  between  the  top  and  bottom  rails  to  maintain  the  correct 
shape.  The  leading  edge  or  simulated  mast  is  formed  by  a  0.750  inch 
diameter  steel  tube.  This  is  equivalent  to  a  mast  radius  of  2.08 
percent  for  the  eighteen  inch  chord.   The  leading  edge  tube  is 
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Figure  3.1  Wind  Tunnel  Model  and  Rig 
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freely  supported  and  removable  once  the  sail  model  is  unrigged. 
This  was  designed  to  allow  for  the  testing  of  several  models  without 
requiring  a  new  rig. 

Mounting  the  rig  and  model  in  the  tunnel  was  accomplished 
by  means  of  a  pin  that  runs  through  the  upper  tunnel  window, 
leading  edge  tube  and  rails,  fitting  into  a  flange  mounted  to  the 
tunnel  floor,  Figures  3.2  and  3.3.  This  method  did  not  require 
extensive  modifications  to  the  tunnel  test  section.  A  truss  type 
structure  downstream  of  the  model  was  incorporated  to  provide  a 
degree  of  torsional  stiffness  to  the  rig.  The  dimension  of  the 
truss  aft  of  the  leading  edge  was  arrived  at  to  coincide  with  the 
breather  slot  downstream  of  the  test  section.  This  allowed  for  the 
linkage  controlling  incidence  to  be  placed  through  the  breather 
slots. 


Figure  3.2   Front  View  of  Model  in  Test  Section 
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Figure  3.3   Side  View  of  Model  in  Test  Section 

C.    WIND  TUNNEL  MODEL  CONSTRUCTION 

The  wind  tunnel  model  that  was  used  in  the  experiment  was 
built  and  donated  by  Trevor  Bayless  at  Waddell  Sails  in  Santa  Cruz. 
The  materials  incorporated  are  the  same  as  those  found  on  a  high 
performance  windsurfing  sail.  The  leading  edge  luff  segment  was 
made  of  a  durable  dacron  cloth  while  the  aft  section  was  made  of  a 
high  modulus  7-mil  mylar  film  sail  material.  A  false  seam  was 
added  to  the  leading  edge  to  accurately  pattern  a  seam  that  is 
present  on  actual  sails.  Areas  beyond  the  body  of  the  sail  section 
were  made  with  nylon  strap  material. 

To  build  camber  into  a  finished  sail  sailmakers  use  a 
combination  of  various  panel  layouts  and  full  length  battens  held 
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in  compression  to  form  the  final  shape.  For  the  tunnel  model  only 
battens  could  be  used  to  form  the  section  shape  since  the  model  is 
2-dimensional.  To  force  the  airfoil  shape  into  the  desired  shape 
preformed  battens  were  used.  The  first  set  of  battens,  formed  from 
0.125  by  0.500  inch  steel,  were  difficult  to  bend  into  the  proper 
shape  and  were  found  to  be  too  flexible  for  use.  This  difficulty 
was  overcome  with  the  use  of  carbon  fiber  -  epoxy  materials. 

To  build  the  carbon  fiber  battens,  a  female  mold  was  first 
made  of  wood  matching  the  desired  camber  less  half  the  estimated 
thickness  of  the  batten.  The  lamination  used  for  the  battens 
consisted  of  the  following  layers  in  an  epoxy  matrix: 

1  6  ounce  S-glass  fiberglass 

2  Carbon  fiber  reinforcing  tape 

1  0.125  inch  balsa  wood 

2  Carbon  fiber  reinforcing  tape 
1    6  ounce  S-glass  fiberglass 

The  final  battens  proved  to  be  very  stiff  and  well  worth  the 
additional  work. 

D.    FLOW  VISUALIZATION 

The  primary  means  of  observing  the  flow  field  around  the 
section  consisted  of  a  series  of  tufts  and  a  smoke  stream.  A 
number  of  different  tuft  materials  and  configurations  were  tried. 
Tape  and  cloth  were  placed  over  the  lacing,  used  to  tension  the 
model,  to  cover  the  ends  of  the  model  to  the  rails.  Flow  along  the 
center  section  was  judged  to  be  adequate  for  the  experiments  with 
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this  configuration.  The  best  tufts  were  determined  to  be  very  fine 
black  thread  with  back  lighting  for  photographs  through  the  side 
windows. 

The  smoke  flow  was  best  observed  through  the  upper  window  with 
flood  lighting  from  the  side  windows.  Reference  8  may  be  consulted 
for  additional  information  concerning  flow  visualization.  The 
smoke  wand  used  in  the  32  by  45  inch  tunnel  was  lengthened  by 
approximately  two  feet  to  minimize  the  growth  of  the  flow  prior  to 
reaching  the  test  section.  It  was  noted  that  the  smoke  flow 
behavior  was  sensitive  to  the  tunnel  velocity.  Velocities  were 
limited  to  the  range  of  30  feet  per  second  for  this  portion  of  the 
experiments. 
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IV.   RESULTS 
A.   COMPARISON  OP  THE  COMPUTED  PRESSURE  DISTRIBUTIONS 

The  first  area  examined  was  comparison  of  the  computed 
pressure  coefficient  (Cp)  distributions  computed  by  the  various 
methods.  The  panel  method  was  found  to  be  particularly  sensitive 
to  the  thickness  of  the  aft  sections.  It  was  desirable  to  keep  the 
thickness  value  as  small  as  possible  to  model  the  actual  physical 
sections.  However,  the  minimum  feasible  value  was  found  to  be  half 
a  percent  of  the  chord  length.  The  Navier- Stokes  solution  would 
have  worked  for  any  thickness  but  a  value  of  0.1  percent  was  used 
to  be  of  the  same  order  as  in  the  panel  method.  Figure  4.1  is  a 
pressure  coefficient  vs  x/c  plot  from  the  panel  code  for  the  11.5 
%  camber  section  at  ten  degrees  incidence. 
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Figure  4.1   Panel  Method  Pressure  Plot 
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The  abrupt  changes  in  the  Cp  aft  of  the  immediate  leading  edge  are 
caused  by  the  change  in  curvature,  from  the  flat  section  at  the 
luff  to  the  second  order  polynomial  curve  at  8.5  %  and  at  the 
transition  to  the  circular  arc  section  occurring  at  26.5  %.  The 
poor  convergence  of  the  pressure  at  the  trailing  edge  is  due  to  the 
wedge  shape  trailing  edge. 

A  comparison  of  the  computed  pressure  using  the  panel,  Euler 
and  Navier-Stokes  methods  is  depicted  in  Figure  4.2  for  ten  degrees 
incidence.  The  difference  in  the  predicted  Cp  approximately  one 
third  aft  of  the  leading  edge  on  the  lower  surface  using  the  Euler 
code  is  caused  by  the  recirculating  flow  region.  A  vector  plot  of 
velocities  for  this  region  predicted  by  the  Euler  calculations  is 
depicted  in  Figure  4.3  while  the  Navier-Stokes  vector  plot  is  shown 
in  Figure  4.4. 
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Figure  4.2   Pressure  Plot  of  Panel,  NS  and  Euler  Methods 
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Figure  4.3   Euler  Velocity  Plot 
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Figure  4.4   Navier-Stokes  Velocity  Plot 
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As  the  incidence  angle  is  increased  the  differences  between 
the  inviscid  panel  code  solutions  and  the  viscous  Navier-Stokes 
predictions  can  be  discerned.  It  can  be  seen  from  Figures  4.5 
though  4.7  that  near  stall  the  panel  code  is  no  longer  valid. 
Figures  4.8  through  4.10  show  the  pressure  field  contours  around 
the  section  computed  by  the  Navier-Stokes  code.  The  high  number  of 
contours  around  the  mast  radius  demonstrates  the  importance  of  the 
large  gradients  and  their  effect  on  the  stall  behavior. 


Pressure  Coefficient 


2.5 1 


a, 
U 


10.0 


Method 
Panel 


Navier-Stokes 


0.0 


02~ 


0.4  0.6 

x/chord 


0.8 


1.0 


Figure  4.5   Panel  and  Navier-Stokes  Method  -  12  degrees 
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Figure   4.6   Panel  and  Navier-Stokes  Plot  -  14  degrees 
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Figure  4.7   Panel  and  Navier-Stokes  Plot  -  16  degrees 
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Figure  4.8  Navier-Stokes  Pressure  Contours  -  12  degrees 
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Figure  4.9   Navier-Stokes  Pressure  Contours  -  14  degrees 
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Figure  4.10   Navier- Stokes  Pressure  Contours  -  16  degrees 
B.   VELOCITY  VECTOR  PLOTS 

Several  significant  phenomena  can  be  observed  by  examining  the 
velocity  vectors  plots.  At  angles  of  attack  below  eight  degrees  a 
well  defined  separation  and  reversed  flow  region  was  predicted  in 
the  area  near  the  lower  luff  with  the  Navier- Stokes  code.  The 
Euler  code  also  predicted  recirculating  flow  but  estimated  the 
region  to  be  significantly  larger  and  to  be  located  further  aft. 
The  reversed  flow  region  was  first  observed  in  the  wind  tunnel  by 
means  of  the  tuft  survey.  The  size  of  the  reversed  flow  region 
predicted  by  NS  was  verified  by  smoke  flow  at  an  angle  of  attack  of 
four  degrees.  Figures  4.11  through  4.13  show  this  behavior  by 
means  of  the  NS  prediction,  smoke  flow,  and  tufts,  respectively. 
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Figure  4.11  Navier-Stokes  Prediction  -  4  degrees 
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Figure  4.12   Smoke  Flow  -  4  degrees 
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Figure  4.13   Tuft  Behavior  -  4  degrees 

At  angles  of  attack  approaching  stall  a  small  separated  flow 
region  was  predicted  on  the  upper  surface  around  the  mast.  This 
region  is  noteworthy  because  the  flow  quickly  reattached  to  the 
flat  luff  segment  and  did  not  cause  the  entire  upper  surface  to 
separate.  This  separation  bubble  was  not  visible  during  the  flow 
visualization  experiment,  nor  was  this  surprising.  Viewing  the 
bubble  was  highly  unlikely  due  to  its  very  small  size.  This 
prediction  is  pictured  in  Figure  4.13. 

A  major  goal  in  the  flow  visualization  effect  was  to  precisely 
identify  the  stall  angle.  It  was  observed,  however,  that  the 
section  did  not  stall  abruptly  in  a  manner  characteristic  of  thin 
airfoil  sections.  The  flow  near  the  surface  appeared  to  slowly 
become  more  turbulent  until  reversing  around  fifteen  degrees. 
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Figure  4.14   Reversed  Flow  Region  -  14  degrees 

This  is  an  approximate  angle  due  to  the  highly  turbulent  and  random 
nature  of  the  flow.  The  Navier-Stokes  solutions  show  the  boundary 
layer  growing  in  thickness  as  the  incidence  is  increased  through 
fourteen  degrees.   This  boundary  layer  growth  is  shown  in  Figures 

4.15,  4. 17, and  4.19  with  accompanying  tuft  photographs  in  Figures 

4.16,  4.18,  and  4.20. 

C.   UNSTEADY  MOTION 

To  increase  the  lift  of  sails  windsurfers  frequently  use  an 
oscillatory  motion  to  increase  the  maximum  lift  generated  by  the 
sail.  An  effort  to  understand  this  phenomena  was  made  using  the 
unsteady  options  in  /ns2.f/.  Several  simulations  were  made  with 
the  test  airfoil  but  only  one  solution  will  be  presented.  To  simulate 
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Figure  4.15   Navier- Stokes  Velocity  Vectors  -  12  degrees 


Figure  4.16   Upper  surface  -  12  degrees 
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Figure  4.17   Navier- Stokes  Velocity  Vectors  -  14  degrees 


Figure  4.18   Upper  surface  -  14  degrees 
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Figure  4.19   Navier-Stokes  Velocity  Vectors  -  16  degrees 


Figure  4.20   Upper  surface  -  16  degrees 
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the  sailors  actions  a  sinusoidal  motion  of  four  degree  about  a  mean 
incidence  of  fourteen  degrees  was  used.  The  reduced  frequency  was 
estimated  to  be  0.19  for  an  oscillation  with  a  1.5  second  period. 
The  code  is  started  from  a  steady  solution  for  time  equal  zero  with 
a  frequency  shift  of  minus  one  half  pi.  The  motion  simulated  is 
shown  in  Figure  4.21  as  a  function  of  time.  The  velocity  field 
that  results  from  the  motion  is  shown  in  Figure  4.22  through  4.26. 
The  velocity  vectors  for  the  simulated  motion  show  attached 
flow  beyond  the  steady  motion  stall  angle.  The  flow  in  fact  stays 
attached  through  eighteen  degrees  but  separates  after  the  motion 
starts  down. 
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Figure  4.21  Unsteady  Motion  Simulation 
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Figure  4.22   Unsteady  Motion  -  11.16  degrees 
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Figure  4.23   Unsteady  Motion  -  13.9  8  degrees 
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Figure  4.24   Unsteady  Motion  -  16.82  degrees 
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Figure  4.25   Unsteady  Motion  -  18.00  degrees 
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Figure  4.26   Unsteady  Motion  -  16.82  degrees 
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V.   CONCLUSIONS 

In  the  computational  and  experimental  investigation  of  the 
flow  field  about  windsurfing  sail  sections  several  important 
findings  were  made.  The  first  area  of  surprise  was  the  extent  of 
flow  separation  that  is  present  over  several  areas  on  the  sections. 
This  was  evident  in  both  the  Navier-Stokes  calculations  and  in  the 
flow  visualization  experiments.  Moreover,  the  degree  of  separation 
was  the  primary  reason  for  the  failure  of  the  two  boundary  layer 
methods  that  were  applied  to  this  problem.  While  the  panel  method 
represents  the  most  elementary  code  used  on  the  problem  it  compared 
favorably  with  the  more  costly  inviscid  Euler  routine.  Since  the 
Euler  code  is  nearly  as  expensive  in  computer  time  as  the  Navier- 
Stokes  method  its  use  can  not  be  justified.  Instead,  much  useful 
information  can  be  obtained  from  the  Navier-Stokes  solutions. 

The  Navier-Stokes  code,  /ns2.f/,  while  performing  well,  has 
been  improved  upon.  Specifically,  a  variable  time  stepping  routine 
has  been  added  which  offers  a  reduction  in  the  number  of  iterations 
and  computer  time  by  over  one  half.  Only  one  turbulence  model, 
Baldwin-Lomax,  was  incorporated  in  the  version  that  was  evaluated 
in  the  thesis.  Several  other  models  have  since  been  added  and 
could  be  evaluated.  Further  investigations  of  the  effect  of  sail 
pumping  used  in  competitive  board  sailing  would  be  of  high  interest 
to  sailors  and  to  aerodynamicists.  Unfortunately,  a  failure  of  the 
computer  system  in  the  CFD  laboratory  prevented  the  modeling  of  a 
realistic  ramp  motion  for  use  in  'pumping'  simulations. 
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The  experiments  undertaken  achieved  their  primary  objective  to 
locate  the  separation  regions.  However,  the  visualization 
experiments  could  be  greatly  improved  upon.  This  could  be 
accomplished  with  the  implementation  of  the  laser  sheet  technique. 
The  rig  and  model  could  also  be  modified  to  include  the  examination 
of  unsteady  motion  effects. 

Finally,  it  is  clear  that  advanced  CFD  methods,  such  as 
Navier-Stokes  solvers,  can  be  used  in  this  field.  Sail  design  is 
a  field  that  thrives  on  new  technology.  The  advances  in  computer 
technology  make  it  only  a  matter  of  time  when  CFD  techniques  will 
migrate  to  the  field  of  sail  design. 


47 


LIST  OF  REFERENCES 


1.  Anderson,  J.  D.,  Fundamentals  of  Aerodynamics.  @d  ed. , 

McGraw-Hill,  Inc.,  1991. 

2.  Cricelli,  A.  S.,  Ekaterinaris,  J.  A.,  Platzer,  M.  F., 

"Unsteady  Airfoil  Solutions  on  Moving  Zonal  Grids, " 
AIAA-92-0543,  Presented  at  the  30th  Aerospace  Sciences 
Meeting  &  Exhibit,  Reno,  NV,  January  6-9,  1992. 

3.  Marchai,  C.  A.,  Aero- Hydrodynamics  of  Sailing.  Dodd,  Mead  & 

Company,  1979. 

4.  Marchai,  C.  A.,  Sailing  Theory  and  Practice.  Dodd, Mead  & 

Company,  19  82. 

5.  Merkle,  C.  L.,  Computational  Fluid  Dynamics  of  Inviscid 

and  High  Reynolds  Number  Flows.  Procopy,  Inc.,  1990. 

6.  Nowak,  L.  M. ,  Computational  Investigations  of  a  NACA  0012 

in  Low  Reynolds  Number  Flows.  Engineer's  Thesis,  Naval 
Postgraduate  School,  Monterey,  CA,  September  1992. 

7.  Smith,  R.  W. ,  "An  Inviscid  Analysis  of  the  Flow  About 

Windsurfing  Sails, "  Presented  at  the  17th  Annual  AIAA 
Symposium  on  Sailing,  Stanford,  CA,  Saturday  and 
Sunday,  October  31,  November  1,  19  87, 

8.  Sommers,  J.  D.,  An  Experimental  Investigation  of  Support 

Strut  Interference  on  a  Three- Percent  Fighter  Model 
at  High  Angles  of  Attack.  Master  of  Science  Thesis, 
Naval  Postgraduate  School,  Monterey,  CA,  September 
1989 


48 


APPENDIX  A:  COMPUTER  PROGRAMS 
Some  of  the  computer  programs  that  were  used  during  the 
research  are  presented  in  the  following  section.  This  listing  is 
by  no  means  a  complete  listing  of  the  software  utilized.  Not 
listed  but  of  significant  utility  are  'hypgen'  and  'plot3d'  both  of 
which  are  well  documented. 
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*****   shape. f 

*  Lt  Matthew  Avila 

*  Thesis  Advisor:  Prof  Platzer 

*  Sail  section  generation  subroutine 

*  Parabolic/circular  arc  or  Lagrange  Interpolation 

*  to  define  camber  line 

*  Wake  points  added  for  hypgen  use       (  hypgen.dat  ) 

*  X-Y  data  points  (  shape. out  ) 
***** 

DIMENSION  X(0:150) ,Y(0:150) ,XNODE(30) ,YNODE(30) 

DIMENSION  xx(300) ,yy(300) ,zz(300) ,Xwake(30) ,Ywake(30) 

OPEN  (101, FILE-' Shape. out' , STATUS =' UNKNOWN' ) 

OPEN  (102, FILE-' hypgen.dat' , F0RM=' formatted' ) 

pi=3. 14159265 

X(0)=1.0 

Y(0)=0.0 

X(150)=1.0 

Y(150)=0.0 

PRINT  *,'This  program  generates  X-Y  section  coordinates  and' 

PRINT  *,'Plot3d  file  for  with  the  HYPGEN  grid  generation 
program. ' 

PRINT  *,'151  points  are  created  for  the  X-Y  plot  and  a 
211x1x1' 

PRINT  *,'grid  in  plot3d  format.' 

10    PRINT  *, 'Enter  camber  line  generation  method' 

PRINT  *,'       1  =  parabola  forward  section  /  circular  arc 

aft' 

PRINT  *,'       2  =  Lagrange  polynomial  defined  by  user' 
READ  *,TYPE 

IF  (TYPE.EQ.l)  THEN 
GOTO  20 
ELSEIF  (TYPE.EQ.2)  THEN 

GOTO  50 
ELSE 

PRINT  *,'You  must  enter  1  or  2 ' 
GOTO  10 
ENDIF 

*  Input  Variables  for  Parabolic-arc  method 

20    PRINT  *, 'Enter  mast  radius  (.02  normal)' 
READ  *,r 

PRINT  *, 'Enter  Maximum  camber' 
READ  *,Ymax 

PRINT  *, 'Enter  Maximum  camber  location' 
READ  * , Xmax 

PRINT  *, 'Enter  Luff  Pocket  Length' 
READ  *,Xluff 

PRINT  *, 'Enter  thickness' 
READ  * , t 
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t  =  t/2 

nfwd=INT(Xluff*70) 

naft=70-nfwd 

rl=( (1-Xmax) **2+Ymax**2) /2/Ymax 

*  Aft  sections  /  Parabolic-Arc 

dxaf t= (1-Xluf f ) /naf t 
DO  30  1=1, naf t 
Xl=I*dxaft 

IF  (Xl.LE. (1-Xmax) )  THEN 
Si=ASIN( (l-Xmax-Xl)/rl) 
Yl=rl*COS (Si) +Ymax-rl 
ELSE 

Yl=Ymax-Ymax*  (  (Xl+Xmax-1)  **2)/  (Xmax**2) 
ENDIF 

X(I)=1.0-X1 
Y(I)=Yl-t 
X(150-I)=1.0-X1 
Y(150-I)=Yl+t 
30     CONTINUE 

*  Wake  section  /  parabolic-Arc 

Si  =  ASIN(  (l.-Xmax)  /  (rl-Ymax)  ) 
Xwake (l)=l+dxaft 
Ywake(l)=-dxaft*TAN(Si) 
DO  40  1=2,30 

dxaft=dxaft*1.2 

Xwake ( I ) =Xwake ( I - 1 ) +dxaft 

Ywake(I)=  -  (  Xwake (I)  -1.  )*TAN(Si) 
40     CONTINUE 
GOTO  100 

*  Input  Variables  for  the  Lagrange  method 

50    PRINT  *, 'Enter  number  of  nodes  to  describe  camber  line' 
READ  *,nn 
DO  60  1=1, nn 

PRINT  *, 'Enter  (x,y)  for  node', I 
READ  *,XNODE(I) ,YN0DE(I) 
60     CONTINUE 

PRINT  *, 'Enter  mast  radius  (.02  normal)' 

READ  *,r 

PRINT  *, 'Enter  Luff  Pocket  Length' 

READ  *,Xluff 

PRINT  *, 'Enter  thickness' 

READ  *,t 

t=t/2 

nfwd=INT(Xluff*70) 

naft=70-nfwd 

*  Aft  sections  /  Lagrange  method 

dxaf t= (1-Xluff) /naf t 
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DO  70  I=l,naft 
XI  =  1  -  I*dxaft 
Y1=0 

DO  80  K=l,nn 
Z  =  YNODE(K) 
DO  90  L=l,nn 
IF  (K.NE.L)  Z=Z*(Xl-XNODE(L) )/(XNODE(K) -XNODE(L)) 

90  CONTINUE 

Yl  =  Yl  +  Z 
80  CONTINUE 

X(I)  =  XI 
Y(I)  =  Yl-t 
X(150-I)  =  XI 
Y(150-I)  =  Yl+t 
70     CONTINUE 

*  Wake  section  /  Lagrange 

100    Si  =  ATAN( (Y(l) +Y(149) )/2/dxaft) 
Xwake(l)=l+dxaft 
Ywake(l)=-dxaft*TAN(Si) 
DO  110  1=2,30 
dxaft=dxaft*1.2 
Xwake(I) =Xwake(I-l) +dxaft 
Ywake(I)=  -  (  Xwake(I)  -1.  )*TAN(Si) 
110     CONTINUE 

*  Lower  Luff  Section  /  Common  for  both  methods 

Thel=pi/2-ATAN(Y(naft) / (Xluf f -r) ) 

X(70)=r+r*COS(thel) 

Y(70)=-r*SIN(Thel) 

dx=X(70) -X(naft) 

slope= (Y(70) -Y(naft) ) /dx 

dxfwd=dx/nfwd 

DO  120  1=1, (nfwd-1) 

X(I+naft)=Xluff+dxfwd*I 

Y(I+naft)=Y(naft) +(X(I+naft) -X(naft) ) *slope 
120     CONTINUE 

*  Upper  Luff  Section  /  Common  both  methods 

The2=pi/2-ATAN( (Xluff -r) /Y(150-naft) ) 

X(80)=r-r*SIN(The2) 

Y(80)=r*COS(The2) 

dx=X(150-naft) -X(80) 

slope=(Y(150-naft) -Y(80) ) /dx 

dxf wd=dx/nfwd 

DO  130  1=1, (nfwd-1) 

X(I+80)=X(80) +dxfwd*I 
Y ( 1+80 ) =Y ( 80 )+(X( 1+80) -X(80) ) *slope 
130     CONTINUE 
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*  Leading  Edge  Section 

Theda=1.50*pi-Thel-The2 

dTheda=Theda/lO 

DO  140  1=1,9 

X (70+1) =r+r*COS (dTheda*I+Thel) 
Y(70+I)=-r*SIN(dTheda*I+Thel) 
140     CONTINUE 


*  SHAPE. OUT  ROUTINE 

DO  150  1=0,150 

WRITE  (101,*)X(I) ,Y(I) 
150    CONTINUE 


*  HYPGEN  ROUTINE 

im  =  211 
jm  =  1 
km  =  1 

j=0 

do  160  i=30,l, -1 

xx ( j ) =Xwake ( i ) 
yy(j)=0. 
zz ( j ) =Ywake (i) 
160    continue 

do  170  i=0,150 
ii  =  i+31 
xx (ii)  =  x(i) 
yy(ii)  =  0. 
zz (ii)  =  y (i) 
170    continue 

do  180  i=l,30 
xx(i+181)=Xwake(i) 
yy(i+181)=0. 
zz(i+181)=Ywake(i) 
180    continue 

rewind  102 

write  (102,*)  1 

write  (102,*)  im,jm,km 

write  (102,*)  (xx(i) ,  i=l,211), 

>  (yy(i),  i-1,211), 

>  (zz(i),  i=l,211) 
close  (102) 

END 
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*****   Shape. out  -  Primary  section  -  11.5  %  chamber 

*****   point  X             Y 

0  1.000000  0.000000E+OO 

1  9.817000E-01  3.269569E-03 

2  9.634000E-01  8.880347E-03 

3  9.451000E-01  1.433349E-02 

4  9.268000E-01  1.963011E-02 

5  9.085000E-01  2.477126E-02 

6  8.902000E-01  2.975801E-02 

7  8.719000E-01  3.459128E-02 

8  8.536000E-01  3.927206E-02 

9  8.353000E-01  4.380124E-02 

10  8.170000E-01  4.817970E-02 

11  7.987000E-01  5.240830E-02 

12  7.804000E-01  5.648780E-02 

13  7.621000E-01  6.041901E-02 

14  7.438000E-01  6.420265E-02 

15  7.255000E-01  6.783941E-02 

16  7.072000E-01  7.133000E-02 

17  6.889000E-01  7.467502E-02 

18  6.706000E-01  7.787510E-02 

19  6.523000E-01  8.093083E-02 

20  6.340000E-01  8.384275E-02 

21  6.157000E-01  8.661140E-02 

22  5.974000E-01  8.923726E-02 

23  5.791000E-01  9.172080E-02 

24  5.608000E-01  9.406247E-02 

25  5.425000E-01  9.626268E-02 

26  5.242000E-01  9.832182E-02 

27  5.059000E-01  1.002402E-01 

28  4.876000E-01  1.020183E-01 

29  4.693000E-01  1.036563E-01 

30  4.510000E-01  1.051545E-01 

31  4.327000E-01  1.065132E-01 

32  4.144000E-01  1.077326E-01 

33  3.961000E-01  1.088129E-01 

34  3.778000E-01  1.097544E-01 

35  3.595000E-01  1.105571E-01 

36  3.412000E-01  1.112212E-01 

37  3.229000E-01  1.117469E-01 

38  3.046000E-01  1.121341E-01 

39  2.863000E-01  1.123831E-01 

40  2.680000E-01  1.124938E-01 

41  2.497000E-01  1.122266E-01 

42  2.314000E-01  1.108858E-01 

43  2.131000E-01  1.084272E-01 

44  1.948000E-01  1.048508E-01 

45  1.765000E-01  1.001566E-01 

46  1.582000E-01  9.434451E-02 

47  1.399000E-01  8.741464E-02 

48  1.216000E-01  7.936696E-02 

49  1.033000E-01  7.020146E-02 
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50  8.499998E-02 

51  8.242778E-02 

52  7.985556E-02 

53  7.728335E-02 

54  7.471112E-02 

55  7.213891E-02 

56  6.956669E-02 

57  6.699447E-02 

58  6.442225E-02 

59  6.185003E-02 

60  5.927781E-02 

61  5.670559E-02 

62  5.413337E-02 

63  5.156115E-02 

64  4.898893E-02 

65  4.641671E-02 

66  4.384449E-02 

67  4.127228E-02 

68  3.870006E-02 

69  3.612784E-02 

70  3.355560E-02 

71  2.842064E-02 

72  2.248211E-02 

73  1.630672E-02 

74  1.048377E-02 

75  5.568937E-03 

76  2.031239E-03 

77  2.082685E-04 

78  2.739891E-04 

79  2.222129E-03 

80  5.866778E-03 

81  9.823438E-03 

82  1.378010E-02 

83  1.773676E-02 

84  2.169342E-02 

85  2.565008E-02 

86  2.960674E-02 

87  3.356340E-02 

88  3.752005E-02 

89  4.147672E-02 

90  4.543338E-02 

91  4.939003E-02 

92  5.334669E-02 

93  5.730335E-02 

94  6.126001E-02 

95  6.521668E-02 

96  6.917334E-02 

97  7.313000E-02 

98  7.708665E-02 

99  8.104331E-02 

100  8.499998E-02 

101  1.033000E-01 


5.991813E-02 

5.618700E-02 

5.245582E-02 

4.872465E-02 

4.499348E-02 

4.126231E-02 

3.753114E-02 

3.379996E-02 

3.006879E-02 

2.633762E-02 

2.260645E-02 

1.887528E-02 

1.514410E-02 

1.141293E-02 

7.681764E-03 

3.950591E-03 

2.194186E-04 

-3.511754E-03 

-7.242926E-03 

-1.097409E-02 

-1.470530E-02 

-1.814092E-02 

-1.984538E-02 

-1.965603E-02 

-1.759094E-02 

-1.384718E-02 

-8.782006E-03 

-2.878779E-03 

3.299166E-03 

9.162276E-03 

1.415104E-02 

1.668940E-02 

1.922775E-02 

2.176611E-02 

2.430446E-02 

2.684281E-02 

2.938117E-02 

3.191952E-02 

3.445788E-02 

3.699623E-02 

3.953459E-02 

4.207294E-02 

4.461129E-02 

4.714965E-02 

4.968800E-02 

5.222636E-02 

5.476471E-02 

5.730307E-02 

5.984142E-02 

6.237977E-02 

6.491813E-02 

7.520145E-02 
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102 

1 

.216000E- 

■01 

8 

.436695E- 

-02 

103 

1 

.399000E- 

•01 

9 

.241464E- 

•02 

104 

1 

.582000E- 

■01 

9 

.934451E- 

-02 

105 

1 

.765000E- 

■01 

1 

.051566E- 

•01 

106 

1 

.948000E- 

■01 

1 

.098508E- 

■01 

107 

2 

.131000E- 

■01 

1 

.134272E- 

•01 

108 

2 

.314000E- 

•01 

1 

.158858E- 

-01 

109 

2 

.497000E- 

■01 

1 

.172266E- 

•01 

110 

2 

.680000E- 

■01 

1 

.174938E- 

•01 

111 

2 

.863000E- 

■01 

1 

.173831E- 

■01 

112 

3 

.046000E- 

■01 

1 

.171341E- 

-01 

113 

3 

.229000E- 

•01 

1 

.167469E- 

-01 

114 

3 

.412000E- 

•01 

1 

.162212E- 

-01 

115 

3 

.595000E- 

•01 

1 

.155571E- 

-01 

116 

3 

.778000E- 

■01 

1 

.147543E- 

-01 

117 

3 

.961000E- 

•01 

1, 

.138129E- 

-01 

118 

4 

.144000E- 

■01 

1 

.127326E- 

■01 

119 

4 

.327000E- 

■01 

1, 

.115132E- 

■01 

120 

4 

.510000E- 

■01 

1 

.101545E- 

-01 

121 

4, 

.693000E- 

■01 

1 

.086563E- 

-01 

122 

4, 

.876000E- 

•01 

1 

.070183E- 

■01 

123 

5, 

.059000E- 

•01 

1, 

.052402E- 

■01 

124 

5. 

.242000E- 

•01 

1, 

.033218E- 

-01 

125 

5, 

.425000E- 

•01 

1 

.012627E- 

•01 

126 

5, 

.608000E- 

•01 

9. 

.906247E- 

•02 

127 

5. 

.791000E- 

•01 

9. 

.672080E- 

•02 

128 

5. 

.974000E- 

01 

9, 

.423725E- 

■02 

129 

6, 

.157000E- 

•01 

9, 

.161139E- 

■02 

130 

6. 

.340000E- 

•01 

8, 

.884274E- 

■02 

131 

6. 

.523000E- 

01 

8, 

.593082E- 

•02 

132 

6. 

.706000E- 

01 

8, 

.287510E- 

•02 

133 

6. 

.889000E- 

01 

7. 

.967501E- 

■02 

134 

7, 

.072000E- 

■01 

7. 

.632999E- 

■02 

135 

7, 

.255000E- 

01 

7. 

.283941E- 

02 

136 

7. 

.438000E- 

01 

6. 

.920265E- 

02 

137 

7. 

.621000E- 

01 

6. 

.541901E- 

•02 

138 

7. 

.804000E- 

01 

6, 

.148781E- 

02 

139 

7. 

.987000E- 

01 

5, 

.740830E- 

02 

140 

8. 

.170000E- 

01 

5. 

.317971E- 

02 

141 

8, 

.353000E- 

01 

4, 

.880124E- 

02 

142 

8. 

.536000E- 

01 

4. 

.427206E- 

02 

143 

8. 

.719000E- 

01 

3, 

.959129E- 

02 

144 

8. 

.902000E- 

01 

3. 

.475801E- 

02 

145 

9, 

.085000E- 

01 

2. 

.977126E- 

02 

146 

9. 

.268000E- 

01 

2. 

.463011E- 

02 

147 

9. 

.451000E- 

01 

1. 

.933349E- 

02 

148 

9. 

.634000E- 

01 

1. 

.388035E- 

02 

149 

9. 

.817000E- 

01 

8. 

.269569E- 

03 

150 

1.000000 

0. 

000000E+00 
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*****   shape. out  -  8.5  %  chamber  section 

*****   point  X  Y 

0  1.000000  0.000000E+00 

1  9.830000E-01  2.340039E-03 

2  9.660000E-01  7.030208E-03 

3  9.490000E-01  1.157152E-02 

4  9.320000E-01  1.596496E-02 

5  9.150000E-01  2.021145E-02 

6  8.980000E-01  2.431186E-02 

7  8.810000E-01  2.826707E-02 

8  8.640000E-01  3.207787E-02 

9  8.470000E-01  3.574507E-02 

10  8.300000E-01  3.926940E-02 

11  8.130000E-01  4.265158E-02 

12  7.960000E-01  4.589228E-02 

13  7.790000E-01  4.899215E-02 

14  7.620000E-01  5.195180E-02 

15  7.450000E-01  5.477183E-02 

16  7.280000E-01  5.745278E-02 

17  7.110000E-01  5.999517E-02 

18  6.940000E-01  6.239951E-02 

19  6.770000E-01  6.466625E-02 

20  6.600000E-01  6.679583E-02 

21  6.430000E-01  6.878867E-02 

22  6.260000E-01  7.064513E-02 

23  6.090000E-01  7.236557E-02 

24  5.920000E-01  7.395033E-02 

25  5.750000E-01  7.539970E-02 

26  5.580000E-01  7.671394E-02 

27  5.410000E-01  7.789332E-02 

28  5.240000E-01  7.893805E-02 

29  5.070000E-01  7.984833E-02 

30  4.900000E-01  8.062432E-02 

31  4.730000E-01  8.126617E-02 

32  4.560000E-01  8.177400E-02 

33  4.390000E-01  8.214792E-02 

34  4.220001E-01  8.238797E-02 

35  4.050000E-01  8.249421E-02 

36  3.880000E-01  8.242350E-02 

37  3.710001E-01  8.205322E-02 

38  3.540000E-01  8.137588E-02 

39  3.370000E-01  8.039147E-02 

40  3.200001E-01  7.910001E-02 

41  3.030000E-01  7.750148E-02 

42  2.860000E-01  7.559588E-02 

43  2.690001E-01  7.338323E-02 

44  2.520000E-01  7.086351E-02 

45  2.350000E-01  6.803672E-02 

46  2.180001E-01  6.490289E-02 

47  2.010000E-01  6.146197E-02 

48  1.840000E-01  5.771400E-02 

49  1.670001E-01  5.365898E-02 
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50  1.500000E-01  4.929688E-02 

51  1.438546E-01  4.589699E-02 

52  1.377091E-01  4.249712E-02 

53  1.315637E-01  3.909725E-02 

54  1.254183E-01  3.569737E-02 

55  1.192728E-01  3.229750E-02 

56  1.131274E-01  2.889762E-02 

57  1.069820E-01  2.549775E-02 

58  1.008365E-01  2.209787E-02 

59  9.469112E-02  1.869801E-02 

60  8.854569E-02  1.529813E-02 

61  8.240025E-02  1.189826E-02 

62  7.625482E-02  8.498384E-03 

63  7.010939E-02  5.098510E-03 

64  6.396396E-02  1.698636E-03 

65  5.781852E-02  - 1 . 701238E- 03 

66  5.167309E-02  -5 . 101112E- 03 

67  4.552766E-02  -  8 . 500986E- 03 

68  3.938223E-02  - 1 . 190086E- 02 

69  3.323679E-02  - 1 . 530073E- 02 

70  2.709139E-02  -  1 . 870059E- 02 

71  2.103181E-02  - 1 . 997337E-02 

72  1.487334E-02  - 1 . 933177E- 02 

73  9.206233E-03  - 1 . 683730E- 02 

74  4.573669E-03  -1 . 272903E- 02 

75  1.419661E-03  -  7 . 400743E- 03 

76  4.650643E-05  - 1 . 363119E- 03 

77  5.858173E-04  4.805154E-03 

78  2.985903E-03  1.051287E-02 

79  7.016724E-03  1.521297E-02 

80  1.229194E-02  1.845497E-02 

81  1.917735E-02  2.024706E-02 

82  2.606275E-02  2.203916E-02 

83  3.294816E-02  2.383126E-02 

84  3.983356E-02  2.562335E-02 

85  4.671897E-02  2.741545E-02 

86  5.360437E-02  2.920754E-02 

87  6.048977E-02  3.099964E-02 

88  6.737518E-02  3.279173E-02 

89  7.426059E-02  3.458383E-02 

90  8.114599E-02  3.637592E-02 

91  8.803140E-02  3.816802E-02 

92  9.491680E-02  3.996012E-02 

93  1.018022E-01  4.175221E-02 

94  1.086876E-01  4.354431E-02 

95  1.155730E-01  4.533640E-02 

96  1.224584E-01  4.712850E-02 

97  1.293438E-01  4.892059E-02 

98  1.362292E-01  5.071269E-02 

99  1.431146E-01  5.250479E-02 

100  1.500000E-01  5.429688E-02 

101  1.670001E-01  5.865898E-02 
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102 

1 

.840000E- 

-01 

6 

.271400E- 

-02 

103 

2 

.010000E' 

-01 

6 

.646197E- 

-02 

104 

2 

.180001E 

-01 

6 

.990288E- 

-02 

105 

2 

.350000E- 

-01 

7 

.303672E- 

-02 

106 

2 

.520000E- 

-01 

7 

.586350E- 

-02 

107 

2 

.690001E- 

-01 

7 

.838322E- 

-02 

108 

2 

.860000E- 

-01 

8 

.059587E- 

■02 

109 

3 

.030000E- 

-01 

8 

.250147E- 

-02 

110 

3 

.200001E- 

-01 

8 

.410000E- 

-02 

111 

3 

.370000E- 

-01 

8 

.539147E- 

-02 

112 

3 

.540000E- 

-01 

8 

.637588E- 

•02 

113 

3 

.710001E- 

■01 

8 

.705322E- 

■02 

114 

3 

.880000E- 

-01 

8 

.742350E- 

■02 

115 

4 

.050000E- 

•01 

8 

.749421E- 

■02 

116 

4 

.220001E- 

-01 

8 

.738796E- 

•02 

117 

4 

.390000E- 

•01 

8, 

.714791E- 

■02 

118 

4, 

.560000E- 

•01 

8 

.677400E- 

■02 

119 

4 

.730000E- 

■01 

8, 

.626617E- 

■02 

120 

4, 

.900000E- 

•01 

8. 

.562431E- 

■02 

121 

5, 

.070000E- 

•01 

8. 

.484832E- 

■02 

122 

5, 

.240000E- 

•01 

8. 

.393805E- 

■02 

123 

5. 

.410000E- 

•01 

8, 

.289331E- 

■02 

124 

5, 

.580000E- 

•01 

8, 

.171394E- 

■02 

125 

5. 

.750000E- 

•01 

8, 

.039969E- 

02 

126 

5. 

.920000E- 

01 

7. 

.895032E- 

02 

127 

6. 

.090000E- 

01 

7. 

.736557E- 

■02 

128 

6, 

.260000E- 

•01 

7. 

.564513E- 

02 

129 

6. 

.430000E- 

01 

7. 

.378867E- 

02 

130 

6. 

.600000E- 

01 

7. 

.179583E- 

■02 

131 

6. 

.770000E- 

01 

6, 

.966624E- 

•02 

132 

6. 

.940000E- 

01 

6. 

.739951E- 

■02 

133 

7. 

.110000E- 

01 

6. 

.499517E- 

02 

134 

7, 

.280000E- 

01 

6. 

245278E- 

02 

135 

7. 

.450000E- 

01 

5. 

.977183E- 

02 

136 

7. 

.620000E- 

01 

5. 

.695180E- 

02 

137 

7. 

.790000E- 

01 

5. 

399215E- 

02 

138 

7. 

.960000E- 

01 

5. 

089228E- 

02 

139 

8. 

130000E- 

01 

4. 

765158E- 

02 

140 

8. 

300000E- 

01 

4. 

426941E- 

02 

141 

8. 

470000E- 

01 

4. 

074508E- 

02 

142 

8. 

640000E- 

01 

3. 

707788E- 

02 

143 

8. 

810000E- 

01 

3. 

326707E- 

02 

144 

8. 

980000E- 

01 

2. 

931186E- 

02 

145 

9. 

.150000E- 

01 

2. 

521145E- 

02 

146 

9. 

.320000E- 

01 

2. 

096496E- 

02 

147 

9. 

490000E- 

01 

1. 

657152E- 

02 

148 

9. 

660000E- 

01 

1. 

203021E- 

02 

149 

9. 

830000E- 

01 

7. 

340039E- 

03 

150 

1.000000 

0. 

OOOOOOE-j 

•00 
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*****  panel. f 

*  Lt  Matthew  Avila  Last  modified  19  Oct  92 

*  Thesis  Advisor:  Prof  Platzer 

*  Panel  Method  Program 

*  Airfoil  shape  input  from  Plot3d  2  dimensional  grid  or 

*  X  -  Y  data  file 
***** 

PARAMETER  (ID=300 , JD=80) 

COMMON  /SET1/X(ID) ,Y(ID) ,DIST(ID) ,SINT(ID) , COST (ID) 
COMMON  /SET2/A(ID, ID) , PI, NPANEL, NODES, NTOT 
DIMENSION  XG(ID, JD) ,YG(ID,JD) 
CHARACTER* 80  FILEIN 
INTEGER  TYPE 
1      FORMAT (A) 

PI=3. 1415926585 

10     PRINT  *, 'ENTER  DATA  TYPE' 

PRINT  *,'        1  =  PLOT3D  FILE' 
PRINT  *,'        2  =  X  -  Y  DATA' 
READ  (*,*)TYPE 
IF  (TYPE.EQ.l)  THEN 
GOTO  2  0 
ELSEIF  (TYPE.EQ.2)  THEN 

GOTO  40 
ELSE 

PRINT  * , ' YOU  MUST  ENTER  1  OR  2 ' 
GOTO  10 
ENDIF 

*  READ   PLOT3D  GRID  POINTS  FOR  AIRFOIL 
20     PRINT  *, 'ENTER  2-D  GRID  FILE  NAME' 

READ  (*,1) FILEIN 

OPEN  (UNIT=201#  FILE=FILEIN,STATUS='OLD'  , FORM= ' FORMATTED ' ) 
PRINT  *, 'ENTER  I  FOR  TRAILING  EDGE  LOWER,  UPPER' 
READ  *,ITEL,ITEU 
REWIND  201 

READ  (201,*)IMAX, JMAX 
PRINT  *,IMAX,JMAX 

READ  (201,*) ( (XG(I, J) ,I=1,IMAX) , J=1,JMAX) , 
>      ( (YG(I, J) ,I=1,IMAX) ,J=1, JMAX) 

1  =  0 

DO  30  J=ITEL,ITEU 

1  =  1  +  1 

X(I)=XG(J,1) 

Y(I)=YG(J,1) 
30     CONTINUE 

NPANEL  =  ITEU  -  ITEL 
NODES   -  NPANEL  +  1 
NTOT    =  NODES  +  1 
GOTO  60 
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*  READ  X-Y  DATA  POINTS  FROM  A  FILE 

40     PRINT  *, 'ENTER  X-Y  FILE  FOR  AIRFOIL' 
READ  (*,1)FILEIN 

OPEN  (UNIT=201,FILE=FILEIN,STATUS='OLD' , FORM= ' FORMATTED ' ) 
PRINT  *, 'ENTER  NUMBER  OF  DATA  POINTS' 
READ  *, NODES 
NPANEL  =  NODES  -  1 
NTOT    -  NODES  +  1 

DO  50  1=1, NODES 

READ  (201,*)X(I) ,Y(I) 
50     CONTINUE 

*  ANGLE  OF  ATTACK  FOR  COMPUTATION  ENTERED 

60    PRINT  *, 'Enter  Angle  of  attack  (degrees)' 
READ  *,AOA 
SINA=SIN(AOA*PI/180) 
COSA=COS (AOA*PI/180) 

*  SLOPE  AND  LENGTH  OF  PANELS  COMPUTED 

DO  70  J=l, NPANEL 
DX=X(J+1) -X(J) 
DY=Y(J+1) -Y(J) 
DIST(J)=SQRT(DX**2+DY**2) 
SINT(J)=DY/DIST(J) 
COST (J) =DX/DIST (J) 
70     CONTINUE 

*  EXECUTE  INFLUENCE  COEFFICIENT,  GAUSS  AND  VELOCITY  SUBROUTINES 

PRINT  *, 'PANELS  COMPUTED' 

CALL  COEFF(SINA,COSA) 

PRINT  *, 'COEFFICIENTS  MATRIX  COMPUTED' 

CALL  GAUSS (SINA) 

PRINT  *, 'GAUSS   REDUCTION  COMPLETE' 

CALL  VELOCITY (SINA, COSA) 

PRINT  *, 'PRESSURE  COEFF.  COMPUTED  (panel. out)' 

END 

****************************************************************** 

*  INFLUENCE  COEFFICIENT  SUBROUTINE 

SUBROUTINE  COEFF (SINA, COSA) 

PARAMETER  (ID=300) 

COMMON  /SET1/X(ID) ,Y(ID) ,DIST(ID) ,SINT(ID) , COST (ID) 

COMMON  /SET2/A(ID,ID) , PI, NPANEL, NODES, NTOT 

DO  110  K=l, NODES 
A (NODES, K) =0.0 
110    CONTINUE 

*  MIDPOINT  NORMAL  VELOCITY=0 

*  I'th  PANEL  MIDPOINT 

DO  120  1=1, NPANEL 
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XMID=0.5*(X(I)+X(I  +  1)  ) 
YMID=0.5*(Y(I)+Y(I  +  1)  ) 
A (I, NODES) =0.0 

*  J  'th  PANEL 

DO  130  J=1,NPANEL 
DLN   =0.0 
BETA  =  PI 

*  I  AND  J  THE  SAME 

IF  (J.EQ.I)  GO  TO  140 
DX1   =  XMID-X(J) 
DX2   =  XMID-X(J+1) 
DY1   =  YMID-Y(J) 
DY2   =  YMID-Y(J+1) 
Dl    =  SQRT(DX1**2+DY1**2) 
D2    =  SQRT(DX2**2+DY2**2) 
DLN   =  LOG(D2/Dl) 

BETA=ATAN2 ( DY2  *DX1 - DX2  *DY1 , DX2  *DX1 +DY2  *DY1 ) 
140  SINTIJ=SINT(I) *COST(J) -COST(I) *SINT(J) 

COSTIJ=COST(I) *COST(J) +SINT(I) *SINT(J) 
A ( I , J) = (BETA*COSTIJ+DLN*SINTIJ) /PI/2 
B= (DLN*COSTIJ-BETA*SINTIJ) /PI/2 
A ( I , NODES ) =A ( I , NODES ) +B 

*  KUTTA  CONDITION  APPLIED  TO  FIRST/LAST  PANEL 
IF  ( (I.GT.l) .AND. (I.LT.NPANEL) )  GO  TO  130 

A (NODES, J) =A (NODES, J) -B 
A (NODES , NODES ) =A ( NODES , NODES ) +A ( I , J) 
130       CONTINUE 

A(I,NTOT)=SINT(I) *COSA-COST ( I) *SINA 
120    CONTINUE 

TEMP=- (COST(l) +COST(NPANEL) ) *COSA- (SINT(l) +SINT (NPANEL) ) *SINA 

A (NODES , NTOT) =TEMP 

END 

********************************************************************* 

*  GAUSS  SUBROUTINE  TO  SOLVE  COEFFICIENT  MATRIX 

SUBROUTINE  GAUSS (*) 

PARAMETER  (ID=300) 

COMMON  /SET1/X(ID) ,Y(ID) ,DIST(ID) ,SINT(ID) , COST (ID) 

COMMON  /SET2/A(ID,ID) , PI, NPANEL, NODES, NTOT 

*  SEARCH  FOR  THE  LARGEST  PIVOT 

DO  200  1=2, NTOT 
IM=I-1 
IMAX=IM 

AMAX=ABS  ( A  ( IM,  IM)  ) 
DO  210   J=I, NODES 

IF  (AMAX.GE.ABS(A(J,IM) ) )  GO  TO  210 
IMAX=J 

AMAX=ABS(A(J,IM) ) 
210    CONTINUE 
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*  SWITCH  ROWS  IF  NECESSARY 

IF  (IMAX.EQ.IM)  GO  TO  220 
DO  230  J    «    IM,NTOT 

TEMP     =    A(IM,J) 

A(IM,J)  -    A(IMAX,J) 

A(IMAX, J)=TEMP 
230    CONTINUE 

*  REDUCE  BELOW  THE  PIVOT  DIAGONAL 
220    DO  250  J=I, NODES 

R     =    A(J,IM)/A(IM,IM) 

DO  240  K=IM,NTOT 

A(J,K)    =   A(J,K) -R*A(IM,K) 
240       CONTINUE 
250    CONTINUE 

*  REDUCE  ABOVE  THE  PIVOT  DIAGONAL 

DO  260  J=IM-1,1, -1 

R     =   A(J,IM)/A(IM,IM) 
DO  270  K=IM,NTOT 

A(J,K)    =    A(J,K) -R*A(IM,K) 
270       CONTINUE 
260    CONTINUE 
200    CONTINUE 

*  CHANGE  PIVOT  VALUES  TO  1 . 0 

DO  2  80  1=1, NODES 
A ( I , NTOT) =A ( I , NTOT) /A ( I ,  I ) 
A(I,I)=1.0 
280    CONTINUE 
END 

****************************************************************** 

*  VELOCITY  DISTRIBUTION  SUBROUTINE 

SUBROUTINE  VELOCITY (SINA, COSA) 

PARAMETER  (ID=300) 

COMMON  /SET1/X(ID) ,Y(ID) ,DIST(ID) ,SINT(ID) , COST (ID) 

COMMON  /SET2 /A ( ID , ID ) , PI , NPANEL , NODES , NTOT 

DIMENSION  CP(ID),Q(ID) 

OPEN  (301, FILE=' panel. out' , STATUS =' unknown' ) 

YMULT=20.0 

DO  310  1=1, NPANEL 
310    Q(I)=A(I,NTOT) 

GAMMA=A (NODES , NTOT) 

*  I'th  PANEL  MIDPOINT 

DO  330  1=1, NPANEL 

XMID=0.5*(X(I)+X(I+1) ) 
YMID=0.5*(Y(I)+Y(I+1) ) 
VTAN=COSA*COST(I) +SINA*SINT (I) 

*  J'th  PANEL 

63 


DO  320  J=1,NPANEL 

DLN=0.0 

BETA=PI 

IF  (J.EQ.I)  GO  TO  300 

DX1   =  XMID-X(J) 

DX2   =  XMID-X(J+1) 

DY1   =  YMID-Y(J) 

DY2   -  YMID-Y(J+1) 

Dl    =  SQRT(DX1**2+DY1**2) 

D2    =  SQRT(DX2**2+DY2**2) 

DLN   =  LOG(D2/Dl) 

BETA=ATAN2 ( DY2  *DX1 - DX2  *DY1 , DX2  *DX1 +DY2  *DY1 ) 
300       SINTIJ=SINT(I) *COST(J) -COST(I) *SINT(J) 

COSTIJ=COST(I) *COST(J) +SINT(I) *SINT(J) 

AA= (BETA*COSTIJ+DLN*SINTIJ) /PI/2 

B= (DLN*COSTIJ-BETA*SINTIJ) /PI/2 

VTAN=VTAN-B*Q(J) +GAMMA*AA 
320       CONTINUE 

CP(I)=1.0-VTAN**2 

WRITE  (301,*)XMID,CP(I) 
330    CONTINUE 
END 
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******  rotate. f 

*  Lt  Matthew  Avila 

*  Thesis  Advisor:  Prof  Platzer 

*  Plot3d  grid  Rotation  Program 
***** 

parameter  (id=300,kd=80) 
dimension  X (id,kd) , Z (id,kd) 
character*15  f ilein, f ileout 

*  Read  in  original  grid 

print  *, 'Enter  plot3d  grid  file' 
10    format (A) 

read  (*, 10) f ilein 

open  (20, f ile=f ilein, form=' formatted' ) 
rewind  20 

read  (20,*)  imax, kmax 

read  (20,*)  ( (X (i, k) , i=l, imax) , k=l,  kmax) , 
>  ( (Z (i,k) , i=l, imax) ,k=l, kmax) 


*  Rotate  the  grid  points 

print  *, 'Enter  the  angle  to  rotate  the  grid  by' 

read  (*,*)theta 

theta=theta*3 . 1415/180 

cost=cos (theta) 

sint=-sin (theta) 

do  30  i=l,imax 

do  30  k=l,kmax 

xold=X(i,k) 

zold=Z(i,k) 

X(i,k)=xold*cost 

Z(i,k) =zold*cost 
30    continue 


-  zold*sint 
+  xold*sint 


*  Write  output  file 

print  *, 'Enter  filename  for  rotated  grid' 
read  (*, 10) f ileout 

open  (40, file=f ileout , form=' formatted' ) 
rewind  40 

write  (40,*)  imax, kmax 

write  (40,*)  ((  X(i,k),  i=l,imax),  k=l,kmax), 
>  ((  Z(i,k),  i=l,imax),  k=l,kmax) 

stop 
end 
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*****  ns2.in  -  input  file  for  ns2 . f 


MACH, 

ALFAO , 

ALFA1, 

ALFARE, 

REDFRE, 

REYNOLDS 

0.10, 

14.00, 

4.0, 

0.0, 

0.19, 

0.80 

ED2X, 

ED2Y, 

ED4X, 

ED4Y, 

ED 

0.00, 

0.00, 

0.030, 

0.030, 

2.0 

DT, 

COUR, 

NITER, 

NEWTIT 

0.0002, 

2100. , 

2160, 

1 

RSTRT, 

OSCIL, 

RAMP, 

NPER 

TSHIFT 

true, 

true, 

false, 

1440, 

-0.5 

TIMEAC, 

IMPLBC, 

EXPLBC 

CIRCOR 

1, 

false, 

true, 

false 

vise, 

BLTM, 

JKTM, 

RNGTM 

true, 

true, 

false, 

false 

ITEL, 

ITEU 

31, 

221 

UNSTST, 

NFRAME 

true, 

8, 

***************************************************************** 
*********** 


Mach 

AlfaO    : 
Alfal 
Redf re 
Reynolds 

ED2x 
subsonic, 


transonic  ) 
ED2z      : 
subsonic, 

transonic  ) 
ED4x     : 
subsonic, 

transonic  ) 
ED4z      : 
subsonic, 

transonic  ) 
ED       : 
ISPEC     : 


Dt 

Cour 
Niter 
Newt it 


Free  stream  Mach  number 
Agle  of  attack,  also  mean  angel  of  attack  for  unsteady 
Amplitude  of  Oscillatory  motion 
Reduced  frequency  k  =  omege  *  c  /  2U 
Reynolds  number  Re  =  cU/n 

X-direction  2nd  order  explicit  smoothing  (  e2x  =  0.00 

0.25  <  e2x  <  0.50 

Z-direction  2nd  order  explicit  smoothing  (  e2z  =  0.00 

0.10  <  e2z  <  0.20 

X-direction  4nd  order  explicit  smoothing  (  e4x  =  0.03 

e4x  =0.05 

Z-direction  4nd  order  explicit  smoothing  (  e4z  =  0.03 

e4z  =  0.05 

Scaling  of  Implicit  smoothing 
Spectral  radious  parameter 

Time  step 

Courant  Number  Cu  =  dt  *  L_max 

Number  of  Iteration  in  this  run 

Newton  subiteration  within  each  timestep 
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RSTRT 
OSCIL     : 
RAMP 
NPER 

dt=T/Nper 
TSHIFT   : 

TIMEAC 
Tacc=0 
IMPLBC    : 
EXPLBC    : 

vise 

inviscid) 
TURBL 
JKTM 
RNGTM 


Restart 

Oscillatory  motion  A(t)  =  AO  +  Al  *  sin  (  k  *  M  *  t  ) 
Ramp  motion 
Number  of  time  steps  in  one  period  of  oscillation, 

Time  shift  in  radiats  to  start  oscilation  for  any  a(t) 

Time  accureat  Tacc=l  and  for  Jacobian  Scaled  Dt, 

Implicit  wall  be  Treatment 
Explicit  wall  be  Treatment 

Viscid  or  inviscid  Boundary  Cconditions  (if  false 

Boldwin-Lomax  eddy  viscosity  (if  false  laminar  only) 
Johnson-King  eddy  viscosity  (not  available) 
RNG  eddy  viscosity  (not  available) 


ITEL,  ITEU  :  Lower  and  Upper  trailing  edge  I  locations 

UNSTST  (Set  true  for  use  with  unsteady  motion  starting  from  an 
steady 

state.  This  will  initialize  time  =  0  for  the  unsteady 
motion 

that  starts  from  a  steady  state  restart. 
NFRAME      :      Number  of  solutions  saved  for  one  cycle  of 
oscillation. 

Saved  starting  at  unit  41. 

Read  grid  from  unit  fort. 11  and  the  flow  from  unit  fort. 31 
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